Accepted for publication in the Astrophysical Journal, November 4, 2006 

Preprint typeset using 1^1^^ style cmulatcapj v. 10/09/06 



RADIATION-HYDRODYNAMIC SIMULATIONS OF COLLAPSE AND FRAGMENTATION IN MASSIVE 

PROTOSTELLAR CORES 

Mark R. Krumholz* 

Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544 



Richard I. Klein 

Astronomy Department, University of California, Berkeley, Berkeley, CA 94720, and Lawrence Livermore National Laboratory, P.O. Box 

808, L-23, Livermore, CA 94550 



O 
O 

(N 

> 
O 

in 

(N 
> 
00 

l> 
O 

o 



o 



X 



Christopher F. McKee 

Departments of Physics and Astronomy, University of California, Berkeley, Berkeley, CA 94720 
Accepted for publication in the Astrophysical Journal, November 4, 2006 

ABSTRACT 

We simulate the early stages of the evolution of turbulent, virialized, high-mass protostellar cores, 
with primary attention to how cores fragment, and whether they form a small or large number of 
protostars. Our simulations use the Orion adaptive mesh refinement code to follow the collapse from 
^ 0.1 pc scales to ^ 10 AU scales, for durations that cover the main fragmentation phase, using three- 
dimensional gravito-radiation hydrodynamics. We find that for a wide range of initial conditions 
radiation feedback from accreting protostars inhibits the formation of fragments, so that the vast 
majority of the collapsed mass accretes onto one or a few objects. Most of the fragmentation that 
does occur takes place in massive, self-shielding disks. These are driven to gravitational instability 
by rapid accretion, producing rapid mass and angular momentum transport that allows most of the 
gas to accrete onto the central star rather than forming fragments. In contrast, a control run using 
the same initial conditions but an isothermal equation of state produces much more fragmentation, 
both in and out of the disk. We conclude that massive cores with observed properties are not likely to 
fragment into many stars, so that, at least at high masses, the core mass function probably determines 
the stellar initial mass function. Our results also demonstrate that simulations of massive star forming 
regions that do not include radiative transfer, and instead rely on a barotropic equation of state or 
optically thin heating and cooling curves, are likely to produce misleading results. 
Subject headings: accretion, accretion disks — equation of state — ISM: clouds — methods: numerical 
— radiative transfer — stars: formation 



1. INTRODUCTION 

The previous generation of telescopes revealed a great 
deal about the gas from which massive stellar clus- 
ters form. With them, observers were able to survey 
the dense clumps of thousands of Mq that are likely 
the proge nitors of clusters, using molecular lin e emis- 
sion (e.g. iPlume et all Il997t iShirley et al.ll2003l). ther- 
mal d ust emission (e.g. iCarev et al.ll2000t Mueller et al.l 
20021) ■ or infrared absorption fe.g. iMenten et al.l l2005t 
Rathborne eFall l2005l . l2006t ISimon et al.ll2006[ ). How- 
ever, the large distances to these regions, their high ex- 
tinctions, and the confusion produced by their density 
prevented these observations from directly probing struc- 
tures with masses comparable to individual stars, data 
that have been available since the 1980s for nearby, low- 
mass star-forming regions. In the last few years, mil- 
limeter interferometers and the Spitzer Space Telescope 
have started to change that situation by making avail- 
able information on the structure of massive star form- 
ing regions comparable to that previously available only 
for low-mass regions. These observations have identified 
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a population of high-mass cores that are comparable in 
mass to individual massive stars. They are dense (mean 
densities ~ 10^ H nuclei cm~^, rising strongly towards 
their centers), cold (temperatures « 10—40 K), turbulent 
(linewidths ~ 1 km s"'^), compact (radii ^0.1 pc), and 
round (aspect ratios o f 2 : 1 or less) (iReid k. Wil son 200^ 
Sridharan et al. 2005; B euther et al.ll2005L I200Q : iGara; 
200 5: Pillai et al. 2006). In some cases they show no 
mid-infrared emission or even mid-infrared absorption, 
indicating that they have not yet converted a significant 
fraction of their mass into stars, and are therefore near 
the onset of star formation. 

The idea that these cores might be the progeni- 
tors of individual massive stars is bolstered by two 
pieces of circumstantial evidence. First, t he mass 
functi on of these cores appears to match the iSalpeteil 
([1955') slope of roughly —1.3 in the logarithmic dis- 
tribution observed for high- mass end of the stellar 
initial mass function (IMF. iBeuther fc Schilkd l200l 
iReid fc Wilso^ l2005l l2006allbl ). This extends earlier 
observations showing that in nearby low-mass star- 
forming regio ns the cor e mass function matches the 
IMF as weU (iMotte et al. 1998: Testi & : Sargent 199^; 
iJohnstone et al.ll2001l : lOnishi et aTEood ). Second, cores 
are mass segregated in such a manner that the mass 
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function is the same throughout a protocluster gas 
clump, with the exception that the most massive cores, 
those greater than sev eral M(d in mass, are found 
only near the ce nter ()Elme green fc Krakowskil [2OOII : 
IStanke et alJ 120 06^. Star clus ters exhibit a very similar 
pattern of mass segregati on ([Hillenbrand fc HartmannI 
[l998t iHuff fc Stahier|[200l . and while some of this may 
be dynamically produced, much of it i s likely a result of 
the locations where t he stars formed ([Bonnell fc Daviei 
[19981 ITan et"ar '2006). 

However, a direct mapping from core mass to star mass 
is only possible if massive cores collapse to form indi- 
vidual stars or smal l - mult iple systems, as proposed by 
iMcKee fc TmI (|2002[ . [200I hereafter MT02 and MT03), 
rather than fragmenting into many objects and produc- 
ing a cluster of low-mass st ars. Whether th i s hap pens 
or not is quite uncertain. [Bate fc Bonnelll (|2005f ) ar- 
gue that dense cores are likely to produce small ob- 
jects because the Jeans mass dec reases with densit y 
at fixed te mperature (although see iMartel et al.ll2006f ). 
iDobbs et~al . (2005) simulate the collapse of massive, tur- 
bulent cores and find that they generally fragment into 
as many as 20 objects, depending on initial conditions 
and on th e assum ed gas equation of state. However, 
iKrumhold (|2006bn uses one dimensional analytic calcu- 
lations to show that radiation feedback from accreting 
protostars can substantially inhibit fragmentation even 
at early times, because at the high accretion rates and 
opacities expected in massive cores, accretion luminosity 
can heat gas to hundreds of Kelvin out to distances of 
^ 1000 AU from an accreting protostar. iKrumholz also 
finds th at using an isothermal or ba rotropic equa t ion o f 
state, as [Bate fc Bon^ (|2005D and lDobbs eFal] (|2005D 
do, is likely to produce misleading results on fragmenta- 
tion because it misses this effect. 

While these calculations are suggestive, because they 
are analytic they are necessarily limited in how they 
deal with real, turbulent cores. The best means of set- 
tling the question of how massive cores fragment is di- 
rect numerical simulation, including a treatment of ra- 
diative feedback from embedded protostars. However, 
such simulations have to this point not been reported 
in the literature. Some simulations of massive star for- 
mation with radiation use quiescent initial conditions 
([Yorke fc Bodenheimeij[l999l : lYorke fc Son nhalter''2002f) 
in two dimensions, and are therefore incapable of answer- 
ing questions about the fragmentation of turbulent struc- 
tures. (Those calculations focus on the effects of radia- 
tion pressure, an important effect in the later evolution 
of massive protostars which we do not consider in detail 
in this paper.) In three dimensions, some simulations of 
massive cores use local cooling func tions rather than solv - 
ing the radiative transfer problem ([Banerjee et al.l[2006f ) , 
and are therefore unable to study the effects of feed- 
back. Moreover, the modifications to the cooling function 
the authors use to approximate the behavior of optically 
thick gas are of unknown accuracy. Simulations of star 
formation with feedback and a treatment of radiative 
transfer have bee n limited to low-mass, no n-turbulent 
initial conditions (Whitehouse fc Bate' '2006") . Further- 
more, both Whitehouse fc Bate and Banerjee et al. only 
advance their simulations to the point where the first 
collapsed object forms, and for this reason they are in- 



capable of studying accretion and fragmentation, or the 
effects of radiative feedback on either of these processes. 

Here we report the first three-dimensional gravito- 
radiation hydrodynamic simulations of the collapse and 
fragmentation of turbulent high-mass protostellar cores. 
Our simulations include radiative transfer and the ef- 
fects of feedback from both accretion onto and nuclear 
burning within embedded protostars. We run through 
the main fragmentation phase and follow the accretion 
process to the point where deuterium burning begins in 
the most massive stars, thereby greatly heating the gas 
and strongly suppressing further fragmentation. This en- 
ables us to address the question of how fragmentation of 
massive cores proceeds, and how radiative feedback in- 
fluences it. In § [2| we discuss the methodology for our 
simulations, and in § [3| we present our results. We dis- 
cuss the implications of these results for the mechanism 
of massive star formation and the origin in of the IMF 
in § [U and summarize in § [5l 

2. SIMULATION METHODOLOGY 

2.1. Evolution Equations 

We describe the evolution of a massive protostellar core 
using the equations of gravito-radiation hydrodynamics 
in the thermal radiation flux-limited d iffusion approxi- 
matio n. Written in conservation form, IKrumholz et al.l 
(|2007( ) show that these equations to leading order in v/c 
are 
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Here p, v, e, and P are the density, velocity, non- 
gravitational specific energy (thermal plus kinetic), and 
thermal pressure of the gas, is the gravitational poten- 
tial, E is the radiation energy density, B = caj^Tg / (47r) 
is the Planck function of the gas temperature Tg, np is 
the Planck-mean specific opacity of the gas measured in 
its rest frame, A and R2 are dimensionless numbers de- 
scribing the radiation field whose significance we discuss 
below, Li and Xi are the luminosity and position of the 
ith star, and n is a unit vector anti-parallel to VE. To 
leading order in v /c, these equations match those of other 
fiux-limite d diffusio n radiation-hydrodynamic codes, e.g. 
ZEUS ( Haves et al.l i200fi'l. 

The gas pressure, specific energy, and temperature Tg 
are related by an ideal equation of state 
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where = 2.33m^f is the mean particle mass in a gas of 
molecular hydrogen and helium with the standard cosmic 
abundance, and we approximate 7 = 5/3 since over most 
of the volume the gas is too cool to excite rotational or 
vibrational modes of hydrogen. In practice the choice 
of 7 has almost no effect, because radiative time scales 
are generally shorter than mechanical ones, so the gas 
temperature and therefore the effective equation of state 
is essentially fixed by radiative transfer effects. 

The gravitational potential is determined by Poisson's 
equation, including the contribution from stars, which 
we treat as point masses: 



VU = 47rG' 
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where Mi is the mass of the ith star. 

The dimensionless numbers appearing in radiation- 
related terms are the flux limiter A and the Eddington 
factor i?2. Their purpose is to interpolate between the 
optically thick and optically thin limits. They are defined 

by 

1^ 1. 
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i^2 = A + A2i?^ (9) 

where kr is the Rosseland-mean specific opacity of the 
gas. The flux limiter has the property that A — > 1/3 
in optically thick regions and A kb.pE/\WE\ in op- 
tically thin regions. For optically thick flows this be- 
havior means that the flux in the frame comoving with 
the gas approaches F —c/{3krp)'VE, the correct 
value for diffusion. For optically thin flows it limits to 
F cEn, so that the effective propogation speed of the 
radiation is limited to c. Similarly, for optically thick 
regions i?2 —>- 1/3, which sets the comoving-frame radi- 
ation pressure tensor to the correct isotropic behavior, 
V — > {E/3)I, where I is the identity tensor. For opti- 
cally thin flows i?2 ^ 1, which gives V Enn, the cor- 
rect limiti ng value for free-s t reami ng radiation. We refer 
readers to iKrumholz et al.l ()2007( ) for a detailed treat- 
ment of the relationship between the comoving frame and 
lab frame quantities, and how the values of A and i?2 are 
related to comoving frame quantities. 

Our equations are easy to understand intuitively. The 
term — AVi? in the momentum equation ([2|) simply 
represents the radiation force K^pF/c, neglecting dis- 
tinctions between the comoving and laboratory frames 
which are smaller than order v/c. Similarly, the terms 
—kpp{4:'kB — cE) and Av • V£' in the gas energy equation 
represent radiation absorbed minus radiation emitted 
by the gas, and the work done by the radiation field on 
the gas. In the radiation energy equation ([4]), the second 
term on the left hand side is the divergence of the radia- 
tion flux, i.e. the rate at which radiation diffuses, and the 
terms on the right hand side describe, from left to right, 
radiation emitted by protostars, radiation emitted minus 
radiation absorbed by gas, work done by the gas on the 
radiation field, and advection of radiation enthalpy by 
the gas. 

Note that our eq uations correspond to those of 
IKrumholz et al.l ()2007f ) for the static diffusion case, which 



IKrumholz et al.l show is the relevant limit for massive 
protostellar envelopes, with two differences. First is the 
addition of the ter ms describing grav ity and point sources 
of radiation, which IKrumholz et al.l do not include. Sec- 
ond is a difference in the coefficient of the work term, 
Av • Vi?. Due to this difference, the equations we give 
here are only accurate to leading order in v/c, rather 
than to first order. In practice, this should make little 
difference in the outcome, since w/c <C 1 everywhere in 
our calculation. 

We discuss the applicability of our equations, including 
the limitations imposed by the approximations we adopt, 
in §0 

2.2. Models for Dust and Protostars 

To complete our specification of the problem, we must 
adopt models to describe the dust, the primary source 
of opacity, and protostellar evolution. We approximate 
that the dust and gas are well-coupled, and neglect the 
possibility that the grain population evolves with time or 
with position except as a function of local gas properties. 
We also assume that dust grains react quickly to changes 
in temperature, so that changes in opacity due to changes 
in the grain population (such as sublimation of certain 
grain components) occur without any time delay. This 
enables us to specify the Planck- and Rosseland-mean 
opacities as simple functions of t he radiat ion tempera- 
ture. We adopt the dust model of iPoUack et all (|l994l ). 
which includes six species of dust grains, each with its 
own sublimation temperature . We approxim ate the tab- 
ulated opacities computed by iPollack et al.l using a sim- 
ple piecewise-linear analytic formula. Our fit gives 
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is in units of K and 
Note that, due to the 



where the radiation temperature T, 
Kp and kr are in units of cm^ g^^ 



optical thickness of a massive core, Tr ~ Tg everywhere 
within one except near its surface. Also note that at high 
temperatures where the dust has sublimed, our choice to 
set Kp = Kr, = 0.1 cm^ g~^ is purely a numerical conve- 
nience we use to represent a "small" opacity. The true 
opacity depends in detail on the radiation spectrum and 
the physical state of the gas (molecular, atomic, or ion- 
ized), but is certainly much smaller than the opacity due 
to dust grains. However, sharp opacity gradients make 
it difficult for our radiation iterative solver to converge, 
so the choice of 0.1 cm^ g~^ is a compromise between 
physical realism and numerical efficiency. This choice 
has little effect in practice, because for the simulations 
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Fig. 1. — Luminosity (upper panel) and radius {lower panel) 
of a protostar versus mass computed using our protostellar model 
with a constant accretion rate of 10"'^ AIq yr~^. The dashed 
line in the upper panel is the luminosity due to accretion. The 
dotted vertical lines mark, from left to right, the masses at which 
deuterium burning starts, deuterium in the core is exhausted, 
convection in the envelope starts, and hydrogen burning starts. 

we describe here only a handful of computational cells 
reach temperatures high enough to be in this regime. 

The final piece of our physical model is a method to 
specify the luminosity of accreting protostars, which ap- 
pears as a source term in the radiation energy equation. 
The input to this model is the mass accretion history of 
the protosta r, which is determined with the sink particle 
algorithm of lKrumholz et ahl (|2004l ) , which we discuss in 
more detail in § 12.31 We adopt the protostellar evolu- 
tion model of MT03, an extension of earlier models by 
[Nakano et al. (1995, 2000). The model is fairly complex, 
so we refer readers to MT03 for a detailed description, 
but we summarize its central features here. The model 
describes a star as a polytropic sphere, and computes 
the evolution of the protostellar radius, central tempera- 
ture, luminosity, and polytropic index using the equation 
of conservation of energy for the star, including terms 
that describe the energy used to dissociate and ionize 
the incoming gas and the energy released by deuterium 
and hydrogen burning. The model includes approximate 
treatments of the onset of deuterium burning in the core, 
the exhaustion of deuterium in the core, the formation of 
a radiative barrier and the formation to a convective en- 
velope, and the start of hydrogen bur ning. It reprod uces 
the detailed numerical simulations of 'Stahlerl (jl988f) and 
[Palla & Stabler (1991) to ~ 10%. Figure [1] shows a sam- 
ple calculation using our numerical implementation of 
the model for the case of a protostar accreting at a con- 
stant rate of 10^^ Mq yr~^. 

2.3. Solution Algorithm 

We solve the evolution equations using the Orion adap- 
tive mesh refinement (AMR) code. The code consists 
of three primary modules, which operate sequentially in 
each update cycle. 

First is the hydrodynamics module 

(jPuckett fc Salt zma iil fTggl iTruelove et all fT998t [KlehU 



Il999f ). which solves the Euler equations of gas dynamics 
([T]), ([2]), and ^ including all of the terms except those 
involving radiation. The hydrodynamics module uses 
a conservative Godunov s cheme with an approximate 
optimized Riemann solver f Torol [l997l ) . The algorithm 
is second-order accurate in time and space for smooth 
flows, and provides a robust treatment of shocks and 
discontinuities using very little artificial viscosity. 

The second part of the code is the gravity module, 
which computes the gravitational potential from the 
Poisson equation ([6t using a multigrid iterati on scheme 
(ITruelove et aLlflOM [KleinllT999l : lFis"heill2002n . 

The third p art is a radiation mod ule, which is up- 
dated using the lKrumholz et all (|2007D operator splitting 
approach, in which we update the dominant radiation 
terms describing diffusion and emis sion minus absorption 
implicitly using the approach of iHowell fc Greenoughl 
([2003), while we update the work and advection terms 
explicitly. This algorithm is stable and accurate for 
problems in the stati c diffusion limit such a s ours , and 
we refer reader s to iHowell fc Greenoughl (|2003D and 
iKrumholz et al.l (|2007D for a detailed discussion of the 
algorithm and the tests we have performed with it. For 
the radiation update, we use the dust model described in 
§ 12.21 to compute the opacities, and we use the luminos- 
ity of each protostar, computed as described below, as a 
source term. 

We supplement these modules by using the Eulerian 
sink particle algorithm of lKrumholz et ahl (|2004f ) to han- 
dle the formation of protostars. When a cell on the finest 
AMR lev el violates the Jeans condition for gravitational 
stability (|Truelove et al.lll997f) . we create a "star parti- 
cle" in that cell. Each such p a rticle is a sink particle, as 
described bv IKrumholz et aLl (j2004D . but also has a cor- 
responding protostellar model, which in addition to the 
mass includes the star's radius, luminosity, polytropic 
index, accretion rate, mass of deuterium remaining, and 
phase of evolution (e.g. whether the star has developed 
a convective envelope yet). In each update cycle, after 
completing the standard update step for sink particles, 
we also update the protostellar model. When we perform 
a radiation update, the protostellar luminosity becomes 
a source term in the radiation energy equation. 

All of these pieces operate with the AMR framework 
(iBerger fc Ohgeilll984l : lBerger fc CollelallTgM iBeiTeTall 
Il994f ). We cover the computational domain with a se- 
ries of levels I — 0,1,2, ... L, where / = is the coars- 
est level, which covers the entire computational domain. 
Each level is a union of rectangular grids, which need 
not be contiguous. The grids are nested, such that ev- 
ery grid on level Z > is entirely enclosed within one 
or more grids on level I — 1. Grids on a given level all 
have the same grid spacing Ax', and spacings on dif- 
ferent levels are related by integer ratios / > 1, so that 
Ax'"*"^ = Ax'' / f. For all the calculations we present here, 
we use f — 2. Each level advances with its own time step 
At', and time steps on adjacent levels obey the relation 
At''^^ — At'//. The process for advancing the calcula- 
tion is recursive. To advance a time step on level I, one 
first updates all the level I cells through a time At' , then 
updates all the cells on level I + 1 through / time steps 
of size At'"'""'^. However, after completing each level I + 1 
update, one advances the cells on level I + 2 through / 
time steps, and so forth down to the finest level present. 
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Fig. 2. — A slice through the xy plane showing the density 
(grayscale) and velocity (arrows) fields at the start of run lOOA. 

After every / cycles on each level Z > 0, we perform a 
synchronization procedure between levels I and I — 1 to 
ensure that mass, momentum, and energy are conserved 
across level boundaries. 

The overall time step is set by the Courant condition 
computed on each level, 



Af = C- 



Aa 



■ Ccff ) 



(12) 



where the maximum is taken over all cells on that level. 
The effective sound speed is 



Ceff 



'7P+ (4/9)£;(l -e-^^P-^^) 
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where 7 is the ratio of specific heats for the gas, and the 
factor (1 — g-i^BP^x^ provides a means of interpolating 
between optically thick cells, where radiation pressure 
contributes to the restoring force and thus increases the 
effective signal speed, and optically thin cells, where ra- 
diation does not provide any pressure. To enforce the 
integer ratio between time steps on different levels, after 
computing the time step Ai' on each level, we set the 
time step on level to At" = min(Ai'/'), then reset the 
time step on all other levels to Ai°//'. 

2.4. Initial, Boundary, and Refinement Conditions 

We chose initial conditions that correspond to the an- 
alytic model of MT03 for high- mass cores. Each calcu- 
lation begins with an initial core that is a sphere of gas 
of mass M and radius ri . The core is centrally concen- 
trated, with density profile p oc r~'^/^, down to some 
inner radius tq. Inside vq the density is constant. This 
corresponds approximately to the thermally-supported 
core and turbulently-supported envelope of a MT03 core. 
To impose the initial turbulent velocity field, we gen- 
erate a 1024^ grid of perturbations with a powe r spec- 
trum P{k) oc k^^ using the method of Dubins ki et al.l 
(1995), corresponding to the spectrum expected for su- 
personic, shock-dominated turbulence. We overlay the 
perturbation cube on the core such that the core just 
fits inside the perturbation cube, and assign the veloc- 
ity at every point in the cube to the corresponding point 
inside the core. The use of a 1024'^ grid allows turbu- 
lent power to be present down to scales of 1/1024 of 
a cloud diameter, which is approximately the diameter 



of the non-turbulent, thermally supported region in the 
MT03 model. We scale the total velocity of the core 
so that the one-dimensional core velocity dispersion is 
a — (GM/ri)^/^, the dispersion required for the core to 
be in approximate hydrostatic balance. We do not drive 
the turbulence or otherwise inject energy after the simu- 
lation starts. We set the initial temperature of the core 
to Tg = 20 K, and the initial radiation energy density to 
1.21 X 10~^ erg cm~'^, the energy density of a blackbody 
radiation field at = 20 K temperature. 

Outside the core is an ambient medium with density 
100 times smaller than the density at the edge of the 
core and temperature 100 times higher, so the core and 
ambient medium are in thermal pressure balance. The 
opacity of the ambient medium is set to zero, so that 
it cannot cool, radiatively heat the core, or inhibit the 
escape of radiation from the core. 

We place the core and ambient medium inside a com- 
putational cube centered on the origin with length L, 
chosen large enough so that no core material ever ap- 
proaches the edge of the computational domain. At the 
boundary we impose symmetry conditions on the hydro- 
dynamic evolution, Dirichlet conditions for the gravia- 
tional field, with the potential on the boundary set equal 
to ~GM/r, and Marshak boundary conditions on the 
radiation field. Marshak conditions are a variant of Neu- 
mann conditions in which the flux into the computational 
domain is set to a constant value of cEo/2, where we set 
Eq equal to the initial radiation energy density. The flux 
out of the computational domain at the face of each cell 
on the boundary is set equal to cE/2, where E is the 
radiation energy density in the cell. This condition im- 
poses a uniform 20 K radiation bath, but allows excess 
radiation generated inside the computational domain to 
escape freely. 

The final piece of our computational setup is the refine- 
ment conditions, which determine when new high resolu- 
tion grids are created or when existing ones are removed, 
a process that occurs automatically throughout the cal- 
culation. We use three refinement criteria. First, in order 
to prevent numerical fragmentation we refine any cell on 
level I in which the density of that cell vi olates the Jeans 
cond ition for self-gravitational stability (|Truelove et al.l 
[1993), 



TTCZ 



G{Ax 



l\2 ■ 



(14) 



We use a Jeans number J = 1/8. Second, we refine any 
cell whose distance from the nearest sink particle is less 
than 16 As', so that the region around sink particles is 
always well-resolved. Third, to ensure that we do not 
produce artificially large radiation pressure forces due 
to poorly resolved radiation energy density gradients, we 
refine any cell in which Ax''\'VE\/E > 0.25, i.e. anywhere 
the radiation energy density changes by more than 25% 
per cell. All of these refinement criteria are applied up to 
a maximum level Lmax, which we specify when we begin 
the calculation. All runs use a resolution of 128"^ cells on 
level and a maxmimum level of 7, giving an effective 
resolution of 16384"^. 

We perform four runs, whose properties we summa- 
rize in Table [T] Run 100 A is our baseline run to which 
we compare the others. We show the initial density and 
velocity field for this run in Figure [2l Run lOOB uses 



6 



Krumholz, Klein, & McKee 



the same parameters as run lOOA and has a turbulent 
velocity field with the same power spectrum, but a dif- 
ferent random realization than run lOOA. This allows us 
to study how the random velocity field influences the re- 
sults. Run 200A uses the same velocity field as lOOA, but 
a more massive core, enabling us to study how our results 
depend on initial core mass. We keep the mean column 
density constant, so that we are not changing the average 
opacity to radiation. Finally, run lOOISO uses the same 
initial conditions as run lOOA, but for it we do not use 
radiative transfer. Instead, we use an isothermal equa- 
tion of state fixed to Tg = 20 K. This amount to setting 
E — and i? = in equations ^ and Q , and dropping 
equation ^ entirely. This lets us isolate how radiative 
transfer affects the results, and study whether simula- 
tions that do not include it are reliable. Run lOOISO is 
very similar i n setu p to some of the models evaluated by 
iDobbs et all ((2001 . 

3. RESULTS 

We evolve each run for 20 kyr, which is 38% of a mean- 
density free-fall time for runs lOOA, lOOB, and lOOISO, 
and 32% of a mean-density free-fall time for run 200A. 
For the radiative runs, in all cases the most massive star 
formed by the end of the run has begun deuterium burn- 
ing, and thus should rise rapidly in luminosity thereafter 
(see Figure [T]), strongly inhibiting further fragmentation. 
The portion of the evolution we follow therefore includes 
the primary fragmentation phase, during which the dens- 
est parts of the core collapse and the pattern of frag- 
mentation is established. As we discuss in § 13.4.21 there 
may be subsequent secondary fragmentation in unstable 
protostellar disks at later times, but this does not signif- 
icantly change where most of the collapsing mass goes. 

In the analysis that follows, we only consider sink par- 
ticles stars if they have a mass of 0.05 Mq or more, the 
mass at w hich "second collapse" to protostellar densi- 
ties o ccurs (jMasunaga et al.iri998l : lMasunaga fc Inutsukal 
[2001 . (Smaller mass objects can still collapse to form 
stars or brown dwarfs, but they produce insufficient pres- 
sure to produce rapid dissociation of molecular hydro- 
gen, leading to second collapse. Instead, they contract 
much more slowly.) We take this precaution because a 
few times over the course of a run our radiation implicit 
solver fails to converge and produces an unrealistically 
low temperature in an isolated cell. The low temperature 
cell may be Jeans un stable and form a sink particle (see 
iKrumholz et al.ll2004 for a discussion of the sink particle 
creation algorithm). These artificially-created sink par- 
ticles are harmless because they contain negligible mass 
- at most a few times 10^^ Mq, usually much smaller. 
Since the radiation solver generally recovers on the next 
time step and produces a reasonable temperature, and 
the regions around the sink particles are gravitationally 
stable once a normal temperature is restored, these low- 
mass sink particles do not accrete or radiate at a not- 
icable rate, and never contain a significant amount of 
mass. We impose a mass threshold before we consider a 
sink particle to be a protostar in order to eliminate these 
spurious objects. This condition may cause us to mis- 
characterize some real, non-numerical but very low mass 
stars, but the amount of mass in stars we could miss in 
this fashion is obviously tiny. Although the isothermal 
run is clearly not subject to this problem, we impose the 



same condition on our analysis of it to avoid introducing 
any bias. 

Each of the radiative runs required roughly 60, 000 — 
70,000 CPU-hours on an IBM SP, running in parallel 
on 128 processors. The isothermal run required approx- 
imately 10, 000 CPU-hours. 

3.1. Summary of Evolution 

We summarize the evolution of our runs starting with 
run lOOA, and we then discuss how the other runs differ. 
We defer discussion of how fragmentation and protostar 
formation occurs across the runs until § 13. 2[ and here 
focus on the overall qualitative morphology of the gas. 

3.1.1. Run lOOA 

Figure [3] shows a time-sequence of the evolution of the 
run, starting from the initial state shown in row (1). Tur- 
bulent motions delay the onset of collapse for a while, but 
as the turbulence decays gas starts to collapse. The first 
object, which we refer to hereafter as the primary star, 
appears 5.3 kyr after the start of the simulation. It forms 
in a shocked filament, which continues to accrete mass 
and by 6 kyr is beginning to form a fiattened protostellar 
disk. Row (2), shows the state of the simulation at this 
point. 

As the evolution continues, several more dense conden- 
sations appear, but most of these are unable to collapse 
and form a protostar before reaching the primary star 
and being sheared apart in the protostellar disk. At 12.2 
kyr a second protostar forms, but it falls into the primary 
star and merges with it at 12.7 kyr, before it has accreted 
0.1 Mq of gas. At the time of merging the primary star is 
already 2.1 Mq, so the mass gained in the merger is neg- 
ligible. Row (3) shows the state of the simulation at 12.5 
kyr, about halfway between when the second star appears 
and when it merges with the primary. We should at this 
point mention a caveat regarding resolution. Due to our 
limited resolution, we are unable to resolve binaries in 
orbits closer than 8 cells, or 60 AU. This means that it is 
likely that at least some of the mergers identified by our 
code are really formations of tight binaries. We discuss 
the implications of this in greater detail in our discussion 
of numerical resolution issues in § 14.21 Here, we simply 
mention that whether the true outcome is a merger of a 
tight binary probably makes very little difference to the 
overall evolution, since the smaller star carries negligible 
mass compared to the primary. 

Only after 14.4 kyr does one of the condensations col- 
lapse to form a second protostar that is not immediately 
accreted, as shown in row (4) . At this point the primary 
star is 3.2 Mq, and has a well-defined massive disk. The 
condenstation from which the new protostar forms is al- 
ready visible in row (3). It is able to collapse and form 
a protostar, unlike several others, because it is fairly dis- 
tant from the central object. This reduces the amount 
of radiative heating to which it is subjected, a topic we 
discuss in more detail in § 13.31 

The next significant change in the system occurs when 
one of the arms of the disk becomes unstable and frag- 
ments to form a third protostar at 17.4 kyr. We show 
the configuration just after this in row (5). At this point 
the central star mass is 4.3 Mq. The disk mass can- 
not be defined precisely, because the disk does not have 



Massive Core Collapse and Fragmentation 



7 



-2, 
0,10 

a.oo 

-Q.05 
-0,10 

o,io| 

0.05 
0.00 
-0.05 
-0,10 
-0,15 

0,10 
0,05 
0.00 
-0.0& 
-0,10 
-0,15 

0,10 
0,05 
0.00 
-0,05 
-0,10 
-0,15 

0,10 
0.05 
0.00 
-0,05 
-0,10 
-0,15 

0,10 
0.05 
0.00 
-0,05 
-0.,1 

-0.15! 



Lag E (g cm ^ 

0-1,0 0.0 1,0 -1-5-0,5 0,5 1-5 -0,5 0.5 1,5 2,5 0,0 1,0 2,0 3,0 



(1) \^ti.Q kyr 



{2) t=&^r 
(4) t=l -^ k yr 



-0,1 O.'T} 0,1 -S -4 + 
pc 1000 AU 




-2 -1 1 -5,0 -2,5 0,0 2,5 
1000 AU 100 AU 



Fig. 3. — Column density as a function of time in run lOOA. From top to bottom, the rows show the cloud state at increasing time, 
as indicated. From left to right, each step to the right corresponds to decreasing the linear size of the region displayed by a factor of 4, 
from a 0.31 pc region in the left column to a 1000 AU region in the right column. At the top of each column we give a scale bar for the 
images in that column. In the left column the region shown is always centered on the origin, and the region shown in the second column 
is indicated by the black box; in the other columns, the region shown is centered on the location of the primary star at that time. Stars 
are indicated by the white plus signs. All images are shown in the same projection. 
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TABLE 1 
Simulation parameters 



Run name 


M (Mq) 


Field 


EOS 


ri (pc) 


ro (AU) 


L (pc) 


Ax-f-."-- (AU) 


pi (lO-l* g cm-3) 


tff (kyr) 


a (km s ^) 


lOOA 


100 


A 


RT 


0.1 


38.4 


0.6 


7.5 


1.0 


52.5 


1.7 


lOOB 


100 


B 


RT 


0.1 


38.4 


0.6 


7.5 


1.0 


52.5 


1.7 


200A 


200 


A 


RT 


0.14 


53.5 


0.85 


10.7 


0.72 


62.4 


2.0 


lOOISO 


100 


A 


ISO 


0.1 


38.4 


0.6 


7.5 


1.0 


52.5 


1.7 



Note. — Col. (3): Perturbation field, A or B. Col. (4): Equation of state, RT = radiative transfer, ISO = isothermal. Col. (7) 
Grid spacing on finest AMR level. Col (8); Initial density of inner, constant density region. Col (9): Free-fall time at the mean 
density. 
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Fig. 4. — Column density at 20 kyr in run lOOA. We do not show the positions of stars. 
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a clearly defined edge. However, we can get an approxi- 
mate mass by defining the disk as all the cells within 1000 
AU of the primary star denser than 10~^^ g cm~'^. This 
gives a mass range of 3.1 Mq. We discuss our definition 
for the disk edge, and disk fragmentation in general, in 
§ 13.41 The fragment is very small compared to the central 
protostar, and remains so as the simulation continues to 
evolve. 

The configuration after 20 kyr of evolution, shown 
in row (6), is substantially similar. Two more small 
disk fragments form, but they both collide with the 
primary star almost immediately after formation, when 
their masses are < 0.1 Mq. This has a negligible effect 
on the mass of the primary. At the end of 20 kyr, the 
primary star is 5.4 M©, the second star is 0.34 M©, and 
the third star, which formed in the disk of the first, is 
0.20 Mq. The disk itself is 3.4 Mq in mass. Thus, the 
system is well on its way to forming a massive star, and 
thus far the vast majority of the collapsed mass has con- 
centrated into a single object. We show a larger plot of 
the full core at this point in Figure S) 

We show the evolution of the primary star's radius, 
mass accretion rate, and luminosity in Figure El The 
model of MT03 gives an accretion rate onto the star-disk 
system of 1.2 x 10"''(m*rf/l M©)^/^ Mq yr~^, where m,^ 
is the mass of the star plus the disk. If we assume that the 
accretion rate onto the star is a fraction m*/m,ci of this, 
we infer an accretion rate of about 2 x lO"** M© yr~^ onto 
the star, which is comparable to the simulation result at 
the end of the calculation; at earlier times, the simu- 
lation gives a higher accretion rate. There are at least 
three reasons for this. First, we assume that the den- 
sity in the central regions was initially constant, whereas 
MT03 assume that the power-law increase in density con- 
tinues all the way to the origin; as a result, accretion is 
delayed for 6 kyr in the simulation, whereas it begins im- 
mediately in the analytic model, thereby spreading the 
accretion over a longer time. Second, MT03 assume that 
the turbulence is undamped, whereas it is damped in 
somewhat less than a crossing time in the simulation, 
which increases the accretion rate. Third, our primary 
star forms out of a single large, shocked filament, and 
accretion onto it is dominated by gas from this filament. 
The shock raises the density and thereby enhances the 
accretion rate relative to what the MT03 model predicts. 

3.1.2. Run lOOB 

Run lOOB uses the same parameters as run lOOA, but 
a different random realization of the initial turbulent ve- 
locity field. The evolution in detail is of course different 
than for run lOOA, but qualitatively it is quite similar. 
Figure [6] shows a time sequence of column densities at 
the same times as in run lOOA. On the largest scales, 
the as in run lOOA, at late times the core is dominated 
by large-scale filaments. There is a primary star at the 
center of this filament network, and on smaller scales it 
is surrounded by a disk, the disk is somewhat denser and 
less extended than in run lOOA, but differs in size by less 
than a factor of 2. 

In run lOOB, a primary star forms at 5.4 kyr, and a 
second object forms at 7.8 kyr. However, the two merge 
at 9 kyr, when the primary is 2.0 Mq and the secondary 
is only 0.17 Mq. By 12.5 kyr, as Figure [5] shows, the 
primary is surrounded by a well-defined flattened disk. 



40 
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Fig. 5. — Radius, accretion rate, and luminosity of the primary 
star versus mass in run lOOA. The dashed vertical line indicates 
the mass at which the star begins burning deuterium. In the 
luminosity plot, the solid line is total luminosity and the dotted 
line is the luminosity due to accretion. The two are identical 
before deuterium burning begins. 

Accretion onto the primary slows down after that point, 
and the disk around it remains fairly stable, until the 
primary collides with a second star at 16.1 kyr. This 
increases the primary's mass from 4.2 Mq to 5.1 M©, 
and the accretion rate increases thereafter, since the col- 
lision brings in an amount of gas considerably larger than 
the amount of mass the primary gains by the collision it- 
self. This also reduces the angular momentum in the disk 
around the primary, allowing more gas to accrete and the 
disk to become more compact. As a result of these two 
effects, the primary grows from 5.1 M© to 8.9 M© in the 
4 kyr after the collision, gaining roughly three times as 
much mass from gas accretion as from the collision. 

The disk around the primary looks similar to that in 
run lOOA at times before the collision at 16.1 kyr, includ- 
ing the beginnings of spiral structure. However, after the 
collision the disk is denser and more compact, and lower 
in mass smaller relative to the star. At 20 kyr, the disk 
mass is 2.4 M©, roughly a quarter of the mass of the star. 
Due to its smaller mass and radius, the disk appears to 
be more stable than the disk in run lOOA, and there is 
no evidence for disk fragmentation. Since this is a result 
of the collision, an event that shows no signs of being 
repeated, it seems likely that further evolution will cause 
the disk to become larger again, and return it to a state 
of instability similar to that of the disk in run lOOA. 

3.1.3. Run 200A 
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Column density as a function of time in run lOOB. See Figure[3]for a detailed descriptic 
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Fig. 7. — Column density as a function of time in run 200A. See Figure |3] for a detailed description. Note tiiat, to accomodate 
tile somewhat larger size of the core, the areas shown in each panel are a factor of 1.5 larger in linear dimension than the analogous panels in[3l 
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In run 200A we use the same turbulent velocity field 
as in run lOOA, but imposed on a 200 Mq core with 
the same initial column density as the 100 Mq core in 
run lOOA. Figure [7] shows the time sequence of column 
densities in run 200A at the same times as in run lOOA. 
The overall appearance of the core is very similar to in 
run 100 A, as are the positions and times of formation of 
the protostars. This is not surprising given the identical 
initial velocity fields. A primary star forms at 5.3 kyr, 
and a second at 11.5 kyr, but this second one merges 
with the primary at 15.1 kyr. At 16.1 kyr, a second 
fragment forms at a larger distance from the primary, and 
it survives. In the final time step, another pair of small 
stars form near it. From 10 kyr of evolution onward, 
the primary has a disk extending several hundred AU 
out. 

At the end of the run, the primary star has reached 
8.6 Mq, and the accretion disk around it is roughly 6 
Mq. As in run lOOA, there is obvious spiral structure in 
the disk. Unlike in run lOOA, there is no disk fragmen- 
tation, although in the final time slice one sees a dense 
condenstation that is analogous to the one that formed 
a protostar out of the disk in run lOOA, and looks likely 
to produce a fragment. The accretion rate onto the star 
is roughly 5 x 10"'* Mq yT~^. This agrees well with the 
MT03 model, which for the initial conditions in run 200A 
predicts an accretion rate of 6.5 x 10^** Mq yr^^ onto a 
14.6 Af© star-disk system. 

3.1.4. Run lOOISO 

We show a time sequence for the evolution of the 
isothermal run at the same times as for run 100 A in Fig- 
ure [8] The overall morphology of the gas is fairly similar 
on large scales, which is not surprising given the identi- 
cal initial conditions. A first condensation forms, and 2 
kyr after the start of the run a primary object forms; it 
remains the most massive star throughout the run. How- 
ever, on smaller scales runs 100 A and lOOISO show very 
significant differences. Due to its lack of thermal sup- 
port compared to run lOOA, the gas in the isothermal 
run is much more filamentary Disks are fiatter, filaments 
have smaller radii, and shock structures are thinner. This 
causes the evolution to proceed quite differently, so that 
by 20 kyr it is fairly difficult to line up the features from 
the two runs except at the grossest level of the location of 
protostellar disks and major filaments. There are many 
more fragmentation sites, and objects are generally much 
clumpier than in the radiative run. The number of pro- 
tostars is considerably larger. 

3.2. Statistics of Fragment Formation 

We summarize the statistics of the stars that form in 
our simulations in Table [2l In addition to describing 
the total number and mass of stars present at the end 
of the run, we give the total number of stars formed, 
the total number of significant stellar mergers, defined 
as those that alter the mass of the more massive merger 
companion by at least 5%, the final mass of the primary 
and of all the other stars combined, and the fraction 
of the primary's mass acquired by mergers rather than 
accretion, defined as the total mass of all stars that merge 
with the primary immediately before the merger, divided 
by the primary's final mass. We show the evolution of 



the mass of the most massive and second most massive 
stars in all the runs in Figure [9l 

The statistics and plot make clear that in all the radia- 
tive runs the primary gains mass almost exclusively by 
accretion. The number of fragments formed is small, and 
the number surviving at the end of the run even smaller. 
Moreover, these fragments always contain an extremely 
small fraction of the total collapsed mass. In none of 
the radiative runs does the primary star gain more than 
~ 10% percent of its mass by collisions. There are kinks 
in the mass versus time curves shown in Figure [9l indi- 
cating sharp rises in the accretion rate, but most of these 
are due to the primary encountering and accreting dense 
gas condensations that had not formed stars, not due to 
mergers. In summary, fragmentation appears to be very 
weak in massive protostellar cores once we take into ac- 
count radiative feedback, and stars appear to gain mass 
by accretion rather than by collisions. 

The fragmentation history is very different in run 
lOOISO. At 20 kyr, there are 7 protostars in run lOOISO, 
and opposed to 3 in run lOOA. There are several more 
that we do not list because they are just below the 0.05 
Mq cutoff. In the isothermal run, these are certainly 
real, since there are no potential problems with an it- 
erative radiation solver. Moreover, a factor of 4 more 
protostars form over the course of run lOOISO than in 
run lOOA, and the fraction of its mass that the primary 
gains by merging rather than accretion is nearly an order 
of magnitude larger. As a result the plot of primary mass 
versus time shown in Figure [9] is much spikier. Some of 
the additional stars form out of the disk around the pri- 
mary, while others form in separate condensations. In 
several cases we can identify analogous condensations in 
the isothermal and radiative runs. In the radiative run 
these are too hot to collapse, and instead they reach the 
primary and are accreted, but in the isothermal run they 
collapse and form protostars. 

The fragmentation we see in run lOOISO is mostly con- 
sistent with previous work on fragmentation in isother- 
mal simulations, which general ly find a great deal of frag- 
ment formation. Dobbs et aQ (f2005l ) find that a simula- 
tion of a 30 Mq turbulent core forms ^ 30 fragments over 
an evolution time slightly longer than ours. They do not 
find any massive stars forming, but that is likely because 
our effective resolution is considerably lower than theirs 
(see i) l4.2p . so they have fewer mergers and more fragment 
formation around their central object. Simulations with 
adaptive SPH codes find that, for isothermal equations 
of state, the amount of fragment formation and the typ- 
ical fragment mass are both highly resolution-dependent 
(jMartel et al.l l2006f l . Runs with radiative transfer are 
unlikely to suffer from this problem, because radiative 
heating shuts off fragmentation on small scales. (It is 
worth noting that we could go to higher resolution in the 
isothermal run, since it is computationally cheaper by a 
factor of ~ 6 — 7 than runs with radiative transfer. How- 
ever, to make the comparison as fair as possible, we use 
the same resolution in the radiative and non-radiative 
runs.) 

3.3. Radiative Heating and Fragmentation 

Clearly there is a significant difference between the ra- 
diative transfer runs, none of which show much frag- 
mentation and for which the morphology is relatively 
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smooth, and the isothermal run, in which numerous frag- 
ments form out of a strongly filamentary morphology 
and a significant fraction of the primary's mass is ac- 
quired through mergers. This suggests that the effective 
equation of state, including radiative heating, is play- 
ing an important role in determining how fragmentation 
occurs. Examining the distribution of temperature and 
Jeans mass, 

in our simulations supports this hypothesis. For refer- 
ence, in the initial state the Jeans mass at the core edge 
is 3.4 Mq, at the mean density it is 2.4 Mq, and at 
the central density it is 0.03 Mq. Note that this Jeans 
mass is 4.71 times the Bonnor-Ebert mass, so centrally- 
condensed objects with masses considerably smaller than 
Mj can still be unstable. 

In Figure [To] we plot for run 100 A the mass M{> T) of 
gas with temperature greater than T, at a time shortly 
after the first fragment other than the primary star forms 
(12.5 kyr) and at the final time in the run (20 kyr), and 
in Figure [TT] we show the spatial distribution of the gas 
as a function of temperature at 12.5 kyr. Clearly by the 
time the second star forms, radiation from the primary 
has heated a significant fraction of the core to well above 
its initial temperature. The temperature is above 50 K 
in 4.4 Mq of gas, including almost all the gas within 
1000 AU of the primary object, which is where much 
of the fragmentation takes place in our isothermal run 
and in other isot hermal simulations in the literature (e.g. 
iBate et al.ll200l . By 20 kyr, the mass heated to more 
than 50 K is 6.0 M©, extending more than 2000 AU away 
from the primary star. 

To study how this is heating is likely to affect frag- 
mentation, in Figure [T2I we show a scatter plot of density 
versus temperature for the mass in run lOOA at 12.5 kyr. 
From the plot, it is clear that almost all of the very dense 
gas, where fragmentation might take place, is heated to 
hundreds of K. However, this heating begins at relatively 
low densities, so almost all the gas denser than 10 times 
the initial density has been heated at least somewhat. 
The p — T distribution is roughly bounded by T oc p"'~^ 
with 7 — 1.2 — 1.3 over the full range of the density distri- 
bution. The exact shape changes at other times, but the 
general feature that the temperature rises continuously 
with density, with no large isothermal density range, per- 
sists at all times after the primary object forms. 

Since the rise is slower than T oc p^/^, the Jeans mass 
does decline as one moves to higher density material, 
and so except at the highest densities and temeperatures 
there is generally more than a thermal Jeans mass of 
material above any given density. At the mean den- 
sity and temperature of the core, there are still many 
thermal Jeans masses. However, the continuous rise of 
T with p in our core still provides an explanati on why 
we s e e so little fragmentation. Simula tions (L i et al.l 
120031 iJappsen et al.jl2005l : iBonnell et all .2006 ) and an- 
alytic work (|Larsonl [20051 ) suggest that how fragmenta- 
tion in a turbulent medium proceeds depends critically 
on the value of 7, with fragmentation proceeding to ar- 
bitrarily small masses as long as 7 < 1, and ceasing for 
7 > 1. The physical argument behind this result is that, 
in a turbulent medium, the densest structures formed by 




Fig. 9. — Mass of most massive star {solid line) and ten times 
tlie mass of the second most massive star {dashed line) as a 
function of time in all runs. Sudden increases in mass correspond 
to points wiiere a smaller star merges with a bigger one. Sudden 
decreases in mass correspond to the points where the title "second 
most massive" star suddenly changes from one star to another 
because the previous second most massive star has merged with 
the most massive. 



the turbulence are generally filaments. Gravitationally 
unstable filaments are able to collapse axisymmetrically 
toward their centers when 7 < 1, but are unable to col- 
lapse axisymmetrically for 7 > 1. If the rate of radiative 
heating and cooling has a density dependence such that 
7 < 1 at low densities and 7 > 1 at high densities, then 
fragments will form at the density corresponding to the 
transition between these two, because at this density con- 
traction of filaments along their axes will stall, and the 
filament will break up into "beads" instead. 

Figure [T^] shows that in a massive core with an ac- 
creting protostar in its center, 7 > 1 effectively over the 
entire core. There is a region of points with 7 « 1, in 
the form of the line of points at T « 30 K at densities 
from 10~^^ — 10~^^ g cm~'^, and indeed these points do 
represent the gas from which the next fragment forms, 
at 14.4 kyr of evolution. They are relatively cool because 
they are ^ 3000 AU from the primary star, and are suf- 
ficiently dense to be self-shielding against its radiation. 
However, these points are the exception, and overall such 
self-shielding distant structures form only rarely. This is 
likely why we see so little fragmentation despite the fact 
that our simulation contains many tens of thermal Jeans 
masses. Filaments do form, but they are unable to even 
begin contracting because radiative heating keeps them 
at 7 > 1. Rather than contracting, stalling, and breaking 
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TABLE 2 
Statistics of Stars Formed 



Run Af20 A'formcd A^mcrgc A^l (Mq) Afothcr (Mq) /merge 

lOOA 3 6 5.4 0.54 0.04 

lOOB 4 7 3 8.9 0.31 0.12 

200A 4 6 2 8.6 0.54 0.06 

lOOISO 7 23 6 7.4 1.5 0.31 



Note. — Col. (2): Number of stars present at 20 kyr. Col. (3): Total 
number of stars formed over the 20 kyr evolution, including those that have 
merged. Col. (4): Number of significant merger events. Col. (5): Mass of 
primary star to 20 kyr. Col. (6): Total mass of all stars but the primary 
at 20 kyr. Col (7): Fraction of primary's mass acquired by mergers. 
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Fig. 10. — Mass of gas above temperature T as a function of 
T in run lOOA. The curves do not reach 100 Mq because this 
analysis only includes gas with a density more than twice the 
initial density at the edge of the core, to ensure that there is no 
confusion with the ambient medium. We show the state of the run 
at 12.5 kyr, just after the second fragment forms {thin solid line), 
and at the end of the run, 20 kyr {thick solid line). We also show 
the distributions at those times computed using temperatures 
derived from a barotropic equation of state rather than the true 
temperatures in our run {thin and thick dashed lines). The top 
axis shows the ratio of Jeans mass Afj at temperature T to Jeans 
mass Mjo in gas of the same density at the initial temperature 
To = 20K. For gas at the mean density of the initial core, 
Mjo = 2.4Mq. 



up into beads, they never begin contracting in the first 
place, and instead their mass drains onto the primary 
star or its disk. 

It is important to note that the energy source in our 
simulation responsible for raising the temperature is al- 
most entirely accretion onto the primary star. At 12.5 
kyr the primary star has not yet started burning deu- 
terium and is only 2 Mq, so the luminosity all comes 
from accretion. At 20 kyr deuterium has ignited but is 
only generating ^ 50% of the total luminosity. Thus, the 
heating does not depend on nuclear burning in the pri- 
mary star. Accretion luminosity by itself is sufficient to 
greatly reduce fragmentation. However, as the luminos- 
ity rises due to nuclear burning, the effect should become 
even more significant. 

For comparison, in Figures fTOl and [TT] we also show re- 
sults using the density distribution in our simulation, but 
using temperatures computed from a barotropic equation 
of state rather than the real temperatures in the simula- 
tion. Simulations have used a variety of barotropi c equa- 
tions of state. We compare to one from (Dobbs et al.l 
l2005f ). which generally produces higher temperatures 



than those used elsewhere (e.g. by iBate et al] 120031 : 
iLi et al.|[200llJappsen et al.li2005D . 

{1, P < Po 

(p/po)'/^ po<p<pi , (16) 
(pl/po)^/^ p>pi 

with To ^ 20 K, po = 10~^'^ g cm^^ ^nd pi = lO^^^ 
g cm~'^. As is clear from the Figures, the barotropic 
equation of state severely underestimates both the tem- 
perature of the gas and the spatial extent of the heated 
region. 

In Figure I12 | we show both t he p ~ T curve result- 
ing from the iDobbs et al.l (|2005f ) barotropic equation of 
state, and also th e curve of the optically thin heating and 
cooling model of iLarsoiil ()2005D . 

^^P>-^'-{ip/por^-\ P>P0' ^^^^ 

with To = 4.4 K, Po = 10^1^ g cm^^, 70 = 0.73, 
and 71 = 1.07. Clearly this equation of state under- 
estimates th e temperature even more severely than the 
iDobbs et al.l bar ot ropic eq uation of state. For either 
the IDobbs et al.l or iLarsonl equations of state, the Jeans 
mass in our simulation is larger than the Jeans mass 
they would predict, often by orders of magnitude. Thus, 
regardless of the details of how fragmentation occurs, 
we expect that simulations that adopt either of these 
proposed equations of state will overpredict the number 
of fragments. Moreover, both the barotropic and opti- 
cally thin equations of state produce a range of density 
with 7 < 1 where the thermodynamics favor fragmenta- 
tion, while our simulation shows that radiation feedback 
largely prevents this type of thermodynamic behavior. 

It is not surprising that the barotropic and optically 
thin cooling eq uations of s tate fa re so poorly. As first 
pointed out bv lKrumhol"3 (|2006bf) . in a collapsing pro- 
tostellar core, before nuclear burning starts the largest 
energy source either internal or external to the core is 
the gravitational potential energy released in the final 
plunge of gas onto the stellar surface. As a result, unless 
one explicitly includes the energy released by accretion 
onto the protostellar surface, and the radiative transfer 
of this energy to the rest of the core, one is ignoring 
the dominant source of energy in the problem. This is 
exactly what the barotropic and optically thin approx- 
imations do. Our results indicate that this is likely to 
result in qualitatively incorrect results for fragmentation 
in massive cores. 

Nor can the problem be fixed simply by using better 
approximate equations of state. As the Figures show, the 
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Fig. 11. — Column density in run lOOA of gas above the temperature indicated in each panel. The top left panel shows all the gas 
(T > 0), and the top row shows gas above temperatures of 50 K, 100 K, and 300 K. The bottom row shows the column density above 
those temperatures, but using temperatures computed from a barotropic equation of state rather than the actual temperature in the run. 
Stars are indicated by white plus signs. The time shown is the same one shown in FigurcO row (3): 12.5 kyr, shortly after the time the 
second protostar forms. 
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Fig. 12. — Scatter plot showing temperature versus density 
for a sample 50,000 cells in run lOOA at time 12.5 kyr. Cells 
are selected with a probability proportional to their mass, so the 
density of points is a true representation of the mass distribution, 
with the exception that we exclude cells with densities below 
twice the initial cloud edge density to ensure that we exclude 
the ambient medium. The diagonal dotted lines are curves of 
constant Jeans mass. The number next to each line indicates the 
value of log(Afj/AfQ) to which it corresponds. The solid lines are 
the curves of p versus T for the barotropic equation of state of 
[Pobbs ct al. (2005) and for the optically thin heating and cooling 
model of .Larson. (,2005i '). 



temperature distribution is a function of both time and 
space, and can change in unexpected ways. For example, 
despite the fact that the luminosity is comparable at 20 
kyr to that at 12.5 kyr, and there is more gas heated to 
moderate temperatures 50 — 100 K, there is actually 
less mass at temperatures ^ 100 K. This is largely be- 
cause an optically thick disk has formed which is shield- 
ing much of the dense gas from protostellar radiation. 
No equation of state that gives the temperature simply 
as a function of the density or other local gas properties 
will reproduce effects like this. 



3.4. Disk Properties 

Here we analyze the properties of the disk around the 
primary star that forms in run lOOA, to better under- 
stand both angular momentum transport and disk frag- 
mentation. To examine the properties of the disk, we 
must first isolate it from the background flow. To do so, 
we choose a density threshold of 10~^^ g cm~'^ to sepa- 
rate disk material from the ambient gas. This threshold 
agrees reasonably well with what one identifies by eye, 
and values up to a factor of ^ 10 different from this do 
not produce qualitatively different results. We also focus 
on gas within 1000 AU of the primary, to ensure that the 
gas we are examining is in orbit around it rather than 
around other protostars. After using these criteria to re- 
move extraneous gas, we compute the total angular mo- 
mentum about the primary of the remaining gas. Since 
our disk is not aligned with the computational grid, to 
analyze it we "deproject" by computing the column den- 
sity E, mass- weighted sound speed c^, and mass- weighted 
angular velocity projected onto the plane orthogonal 
to the angular momentum vector. We do this at 17 kyr, 
just before the first disk fragment forms, and at 20 kyr, 
the end of our run. The deprojected maps of disk prop- 
erties at these times form the basis of our analysis. We 
show the disk at these two times, before deprojection, in 
Figure [H 



3.4.1. Angular Momentum Transport 

First we wish to determine the effective viscosity a for 
our disks. In a Keplerian disk this i s related to the kine- 
matic viscosity uhy u = 2acl / (Sfi) (jShakura &: SunvaevI 
IT973h . The inward radial drift velocity of the material 
in the disk is vn ~ 1^/ R, and the accretion rate onto the 
central star is M = 2t: RY,vr, so a « MQ/ {Snivel). We 
therefore estimate a by computing the mass-weighted 
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Fig. 13. — Column density of the disk in run lOOA at 17 kyr [top row) and 20 kyr {bottom row), in two orthogonal projections. Stars 
are indicated by white plus signs. 
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For the properties of our disk at either 17 kyr or 20 
kyr, and taking M to be the mean accretion rate be- 
tween these times, 4 x 10"'* Mq yr~*, we find an ef- 
fective a ~ 1.0 — 1.6. This is quite rapid angular mo- 
mentum transport, and is significantly larger than what 
one expects due to purely local transport phenomena 
(e.g. lGammiell200ll ). Accretion is obviously highly time- 
dependent and unsteady, and the disk never settles into a 
steady state, so the rate of angular momentum transport 
at any instant may be very different from the average. 

To understand the angular momentum transport 
mechanism, we analyze the spiral pattern in the disk 
by computing Fourier coefficients of the density distri- 
bution in azimuth around the primary star. Defining r 
as the distance from the primary star in our projection 
and (/) as an angular coordinate in the projection plane, 
we compute 



dre™'^rE(r,( 



(19) 



We plot the normalized power |cmP/|co|'^ for m = 1 — 10 
in Figure [Ml As the figure shows, at both 17 kyr and 
20 kyr the vast majority of the power is in the m = 1 
spiral mode, with smaller amounts of power in other odd 
modes and very little power in even modes. This sug- 
gests that angular momentum transport and spiral arm 
formatio n in our disk is primarily due to the SLING in- 
stability (|Adams et al.lll989t IShu et al.lll990[ ). For disks 
as massive as ours, Md ~ 0.5M*, this mechanism enables 
accretion on a disk dynamical time scale rather than a 
viscous time scale, consistent with our value of a ~ 1. 
Angular momentum transport occurs via a global rather 
than a local instability. Our result is also consistent 
with previous si mulations of gravitational instab ility in 
massive disks bv lLaughlin fc Bodenheimerl ()1994[ ). which 
show that most power goes into the to = 1 mode. 

The behavior of disks in our simulations is signifi- 
cantly differen t than tha t seen in t he si mulations of 
iLodato fc Rici (|2005D and lRice et all (|2005D . who model 
the evolution of disks with masses ranging from 10% to 
100% of the primary object mass using polytropic equa- 
tions of state with 7 > 1, with added cooling terms that 
remove energy on time scales from 3 — 13il~^. They 
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Fig. 14. — Normalized power |cmp/|cop in azimuthal mode m 
in the disk around the primary star, at 17 kyr {solid line) and 20 
kyr {dashed line). 

find that for all their runs angular momentum trans- 
port is primarily local and that accretion occurs on a 
viscous rather than a dynamical time scale, with values 
of a < 0.06 in all stable disks. When spiral arms form, 
the majority of the power is in m = 2 modes. Disks settle 
into steady states except when Md ^ . 

This difference in behavior is likely to be a real physical 
eff'ect, caused by two diff erences between the properties 
of our disk and those of iLodato fc Ricd and I Rice et all 
First, the temperature in our disk is set almost entirely 
by radiative heating and cooling, in contrast to the poly- 
tropic plus cooling runs in which disk temperatures are 
set by the balance between viscous heating and radiative 
cooling. Viscous heating does not provide enough energy 
to raise the disk temperature significantly in our run. As 
a result, the temperature in our disk is almost entirely a 
function of distance from the primary star, with no sig- 
nificant variation at a given distance due to spiral arms 
or other density or velocity structures in the disk. Conse- 
quently, the thermodynamic behavior of our disk is closer 
to an isothermal equation of state than to the stiffer equa- 
tions of state produced by values of 7 > 1 and cooling 
time scales longer than the disk orbital period. This fa- 
vors the growth of large-scale global modes that produce 
rapid angular momen tum transport, an effect pointed 
out by ILodato fc Ricd to explain the differences between 
their s imulations and those of [Laughlin & Bodcnheime^ 
11994). 

Second, in our simulation the disk is never stable or iso- 
lated. The average accretion rate of 4 x 10^^ Mq yr~^ 
from 17 to 20 kyr, assuming all mass that reaches the star 
is processed through the disk, corresponds to ~ 30% of 
the disk mass per orbital period. This is obviously a 
huge perturbation. Most of this accretion comes from a 
large filament (as shown for example in Figure [4]) that 
is sheared out into a disk as it approaches the protostar, 
an effect that obviously favors m = 1 spiral structure. 
Partly as a result of this perturbation, our disk is never 
able to settle into a quasi-steady state, and it forms sev- 
eral fragments. These likely aid in shepherding material 
inward into the primary star. 

These two effects suggest t hat our results are no t in- 
consistent mththe findings of lLodato fc Ricd (|2005| ) and 
iRice et al] ()2005D . simply that our disk is in a different 
regime of parameter space than they have explored. 
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Fig. 15. — Deprojected column density {upper panel) and 
Toomre Q {lower panel) of the disk around the primary star in 
run lOOA at 17.4 kyr, just as the first disk fragmentation occurs. 
The positions of stars are indicated by the white plus signs. In 
the plot of Q, the black contour indicates Q = 1, and the exterior 
black region consists of points for which there is no gas above the 
density threshold we use to define the edge of the disk. 

3.4.2. Disk Fragmentation 

To hel p unders t and w hy our disk fragments, we com- 
pute the iToomrd ()1964[ ) parameter 



Q 



ttGE' 



(20) 



We do this at each point, and we also compute the mass- 
weighted azimuthal average 



{Q)Ar) 



/:^d0S(r,0) 



(21) 



as a function of radius. Obviously this calculation is 
somewhat approximate, since we do not have a thin, 
steady, Keplerian disk with a well-defined edge as a back- 
ground state. Nonetheless, it can give us some insight 
into the state of the disk and the reason that it frag- 
ments. 

Figures [15] and [16] show plots of deprojected column 
density and Q for the disk in run 100 A at a time of 17 
kyr, just before the first fragment forms in the disk, and 
at 20 kyr, when the disk-formed star is still present but 
there are no other obvious fragments forming. We show 
the azimuthally-averaged column density and Toomre Q 
as a function of radius at these two times in Figures [T7] 
and [18] As the plots show, at 17 kyr there is a broad 
region where Q < 1, both at individual points and in 
the azimuthal average. This corresponds to the approx- 
imate location where the fragment forms. At 20 kyr, it 
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Fig. 16. — Deprojected column density {upper panel) and 
Toomre Q {lower panel) of the disk arou nd the primary star in 
run lOOA at 20 kyr. For details see Figure [TSl 
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Fig. 17. — Azimuthally-averaged column density {upper panel) 
and Toomre Q parameter {lower panel) as a function of radius in 
the deprojected protostellar disk at 17 kyr, just before the first 
disk fragment forms. 



is clear that there is another region of the disk that has 
Q < 1 at individual points, but there is no region where 
the azimuthally-averaged value of Q < 1 (although the 
ring at 350 AU is extremely close). In general, the disk 
seems more stable than at 17 kyr, which likely explains 
why there is no obvious fragment formation at 20 kyr. 
However, it seems entirely possibly that further evolu- 
tion would decrease Q and lead to additional fragment 
formation. Although the luminosity of the star will con- 
tinue to rise as it gains mass, heating the disk, the disk 
mass will also rise, increasing its shielding against pro- 
tostellar radiation. It is not clear which of these effects 
will dominate. 
While fragment formation in massive protostellar disks 
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Fig. 18.— Same as Figure [TSl but at 20 kyr. 

is an interesting phenomenon, it seems unlikely to be a 
significant hindrance to accretion onto the primary star 
at this point in the evolution. There is no noticable drop 
in the accretion rate onto the primary after 17 kyr, and 
from 17 kyr to 20 kyr the primary star mass increases 
by a factor of 6 more than the mass of the embedded 
fragment. Clearly, most of the mass in the disk is going 
into the primary, not into disk-formed fragments. 

Our findings on disk fragmentatio n are broadly consis- 
tent w ith the analytic predictions of iKr after fc Matzneil 
(|2006( l. who analytically model massive protostellar disks 
and find that they are unstable to fragmentation at 
radii ^150 AU for central stars of mass >4 
iKratter fc Matzneil . extending work of iMatzner fc LevinI 
(|2005l ). predict that steady-state disks should fragment 
if their sound speeds fall below a critical value Ccrit ~ 
1.04(GM)^/'^, where M is the accretion rate onto the 
star-disk system. 

Before applying this condition, we modify it in two 
ways. First, rather than using the isother mal sound 
speed to compare to Ccrit as lMatzner fc LevinI suggest, we 
use the adiabatic sound speed because our disks are opti- 
cally thick to their own radiation. Second, we modify the 
criterion to account for the fact that angular momentum 
transport in our disks appears to be du e to a global rather 
than a local gravitation al instability. IMatzner fc Leviiil 
and IKratter fc Matznerl determine disk stability using 
the criterion of iGammid (pOOl ). who simulates angular 
momentum transport by local gravitational instabilities 
and finds that these produce a maximum effective vis- 
cosity a = 0.23. The critical sound speed depends on a 
as Ccrit oc a^/"^, because a determines the rate at which 
material is processed through the disk onto the central 
object. Consequently, the increased rate of angular mo- 
mentum transport in our disks make s disks more stable. 

With these two modifications, the IKratter fc Matzneil 
critical temperature for instability to fragment formation 
becomes 



;rit«0.41-^(aGM)2/3 



= 39a^/^ 




(22) 
(23) 



For the mean accretion rate of 4 x 10 Mq yr ^ from 17 
kyr to 20 kyr, and our range of estimates a — 1.0 — 1.6, 
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this gives Tcrit w 70 — 100 K. The temperature in the 
outer parts of our disks is generally in this range, which 
explains why they are marginally unstable to fragment 
formation. It also explains why our typical fragmenta- 
tion radius is somewhat larger than iKratter fc Matzneil 
predict. 

4. DISCUSSION 

4.1. Discussion of Physical Approximations 

Our physical formulation of the problem contains three 
significant simplifications. Here we discuss them, with 
the goal of assessing how much they might affect our 
results. 

First, our treatment of radiative transfer, though a sig- 
nificant improvement on previous three-dimensional cal- 
culations that ignored radiation entirely, is still quite ide- 
alized. Our approach is gray, so we miss effects that arise 
from the frequency-de pend ent opacities of dust grains. 
iPreibisch etall (|1995D and lYorke k Sonnhalteil (|2002D . 
based on two dimensional calculations, find that, com- 
pared to gray, multi-frequency calculations generally pro- 
duce higher dust temperatures and greater degrees of 
anisotropy in the radiation field. The fact that, even 
with the anisotropy effect, gray radiation almost always 
underestimates the true temperature suggests that our 
results on fragmentation are fairly secure, since higher 
temperatures would further reduce the level of fragmen- 
tation. Nonetheless, increasing the anisotropy of the ra- 
diation could could conceivably leave some parts of the 
flow cooler than we find, so it would be useful to repeat 
some of the calculations we present here with a multi- 
frequency radiative transfer code. 

In addition to being gray, our radiative transfer ap- 
proach uses the flux-limited diffusion approximation, 
which is an approximation that is only highly accurate 
in regions that are very optically thick (or that are very 
optically thin and for which the geometry is simple) . Our 
cores have initial mean surface densities of 0.66 g cm~^, 
which gives them optical depths 1 to infrared photons. 
Thus, the cores overall are marginally optically thick on 
average. However, in the dense regions where fragmenta- 
tion takes place and where it is most important to treat 
the radiation correctly, surface densities are more typi- 
cally tens or hundreds of g cm~^. Thus, the regions in 
which fragmentation occurs are extremely well-described 
by the diffusion approximation. Thus, we consider it 
extremely unlikely that our results would change quali- 
tatively if we were to use a more accurate treatment of 
the radiation field. 

A second limitation to our approach is the uncertainty 
in what our initial and boundary conditions should be, 
based on our imperfect knowledge of the properties of 
massive cores and their environments. The MT03 model 
which we have used fits the data reasonably well, but ob- 
servations to date reveal only a little about the internal 
structure of massive cores. The primary uncertainties 
likely to affect our results are the degree of central con- 
centration, the amount of internal turbulent structure, 
and the nature of any external perturbations. The first 
of these is difficult to determine because observations of 
massive core gas are generally made with interferome- 
ters such as the SMA, which remove large-scale power 
and thus make it difficult to determine quantities like 



large-scale density gradients. The amount of internal 
turbulent structure is poorly known observationally sim- 
ply because of resolution and sensitivity limits. Massive 
cores are too small for observations to determine fine de- 
tails of their internal structure. The nature of external 
perturbations is uncertain because massive cores are em- 
bedded within clouds that are themselves turbulent, so 
rather than providing a constant pressure boundary, the 
external environment may fiuctuate and drive turbulent 
motions in a massive core. 

One might expect that the stronger the central con- 
centration, the less fragmentation will occur. However, 
non-radiative models produce a great deal of fragmen- 
tation regardless of the initial degree of central concen- 
tration, even for cases much more conc entrated than the 
k„ = 1.5 density gra dient we use fe.g. iBate et al.l[2003t 
iGoodwin et all (2004 iDobbs et al.ll2005D . Thus, the de- 
gree of initial concentration in our models seems unlikely 
to change our results qualitatively. For turbulent den- 
sity structure, whether it is structure present initially or 
structure induced by external perturbations on a core's 
surface, obviously a higher degree of internal structure 
favors more fragment formation, and massive cores are 
likely to have some internal structure because they are 
turbulent. Thus, we do not regard our finding on the 
absolute number of fragments formed to be definitive for 
what will happen in real high-mass cores. However, since 
we use identical initial conditions for the radiative and 
isothermal runs, the differences between them are likely 
to remain even for more structured or less concentrated 
initial density fields. Morever, the general effect we have 
found, that radiative heating can shut off collapse to 
small fragments, while allowing the first object formed 
to grow to large masses by accreting gas that cannot 
otherwise collapse, should apply regardless of the initial 
and boundary conditions. Thus, although the quanti- 
tative results may vary for different initial or boundary 
conditions, the qualitative results that massive cores do 
not fragment very strongly should be robust. 

The third major limitation to our approach is that 
we have neglected magnetic fields. This simplification 
is partly justified by observational ignorance. Even in 
nearby low-mass star-forming regions there is consider- 
able controversy ove r how dynamically significant mag- 
netic fields are (e.g. lCrutchejfl999t iPadoan et al.l[200l 
iTassis fc Mouschoviasll2004D . and observations of more 
distant, obscu red high-mass st ar forming regions are far 
more difficult. iCriitched (l2005l ) reviews the available data 
and concludes that magnetic fields are marginally dy- 
namically significant, but this conclusion is highly un- 
certain due to potential systematic errors in transform- 
ing_the observed signal into a magnetic field strength (see 
lKrumholzii2006a for a discussion of this point). Even if 
magnetic fields are dynamically significant in the initial 
core, simulations show that turbulence ca n significantly 
accel erate the rate of ambipolar diffusion (jHeitsch et al.l 
|2004[ ). so the field might diffuse out quickly and have 
little effect on the overall evolution. Moreover, even if 
magnetic fields are dynamically important and can mod- 
ify how fragmentation proceeds, our result that radiative 
transfer suppresses fragmentation should still hold qual- 
itatively. 

A final magnetic effect that could be significant is on 
our protostellar disks. Since we have no magnetic fields. 
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we obviously have no magnetorotational instability and 
its associated angular momentum transport. It is unclear 
if the MRI can operate in massive protostellar disks, be- 
cause their column densities of 100 g cm~^ may render 
them so opaque to protostellar ultraviolet radiation that 
their ionization fractions will be too low for the MRI to 
operate. (MRI requires sufficient ionization for the mag- 
netic Reynolds n umber to be greater than about unity, 
iSano et aT]|1998l ) However, if the MRI does operate, its 
effect should be to increase the rate of angular momen- 
tum transport, which will raise the mass of the primary 
and lower the surface density of our disks, thereby reduc- 
ing their propensity to fragment. Thus, if anything we 
overestimate fragmentation in our radiative runs. 

4.2. Numerical Resolution 

Another potential concern is whether our results de- 
pend on our numerical resolution. At some level, they 
probably do. The accretion radius around our sink par- 
ticles is 4 cells, which corresponds to 30 AU in runs lOOA, 
lOOB, and lOOISO, and 43 AU in run 200A. This means 
we are unable to resolve binaries whose closest approach 
is smaller than twice this, and stars that should become 
tight binaries will instead be merged in our code. How- 
ever, given that except in the isothermal run, the amount 
of mass gained by mergers is negligible, this effect cannot 
be significant in setting the primary star mass. We are 
also insensitive to the formation of fragments on scales 
smaller than the accretion radius, so we could conceiv- 
ably underprcdict the number of fragments. In our radia- 
tive runs this is unlikely to be a problem, because, as Fig- 
ure [Tl] shows, all the gas within ^ 30 AU of the primary 
protostar is heated to ^ 300 K even at very early times. 
Thus, heating should very strongly suppress fragment 
formation there. On the other hand, this effect could 
be very significant in the isothermal run, since isother- 
mal calculations with higher resolution than we have 
used do find significant amounts of fragment formation 
on ^ 100 AU scales (e.g. iBate et al.l [20031 iDobbs et all 
l2005l ). and simulations with varying resolution find that 
the amount of fragmentation and the m ean fragment 
mass are resolution-dependent (jMartel et al. 2006). 

A related concern is that we might be missing fragmen- 
tation in our disks due to excessive numerical viscosity. 
Our Cartesian grid pr oduces a numerical vis cosity that 
gives an effective a of (jKrumholz et al.ll2004f ) 



a«78 



-3.85 



100 AU 



-3.85 



(24) 
(25) 



where r is the distance from the star about which the 
disk orbits, the mass of the star is Mi in units of Mq, 
the disk temperature in units of 100 K is Tipo, and Aa;7.5 
is the grid spacing in units of 7.5 AU. This means that, 
at distances ~ 300 AU, where much of the disk mass re- 
sides, our typical numerical a is of order 10^^, insignif- 
icant compared to that induced by the SLING instabil- 
ity. On scales ^ 100 AU, however, numerical viscosity is 
probably significant in shaping the evolution of our disks, 
and may inhibit fragment formation. This problem, to 
the extent that it is significant, almost certainly affects 
the isothermal run more than the others, since irradiated 



disks at such small radii are quite hot and thus resistant 
to fragmentation. The problem may also be more severe 
in run 200A, where the larger cell size means that the 
untrustworthy region is a factor of 2 — 3 larger. This 
may explain the reduced disk fragmentation we find in 
that run. 

4.3. Implications for Massive Star Formation 

One of the ouststanding problems in massive star for- 
mation has been how to gather enough mass to make a 
massive star. Our simulations, coupled with other re- 
cent work, suggest that we may be nearing a solution. 
Observations now unambiguously reveal that there are 
massive, cen trally concentrated cores (e.g. see review by 
lGarayl[2005l ) . How these cores form is a topic of active re- 
searchie^lPadoaniLNordlundl Hool; iTillev fc Pudritj 
[2003lLi et al. 2004), but since we can determine massive 
core properties from observation we need not solve this 
problem to model their evolution. 

While it is appealing to see massive cores as the progen- 
itors of massive stars, one might legitimately worry that 
these objects fragment to produce large numbers of low- 
mass stars rather than a few massive stars. Some of these 
low-mass stars c ould possibly gain m as s via competitive 
accretion of gas (Bonncll et al. 2001a, b, 2004) or via stel- 
lar collisions (BonneU et al. 1998; Bonncll & Bate 2005), 
becoming massive stars. However, the former mechanism 
seems not to work unless star cl usters are typically sub- 
virial, globally collaps i ng objects ( Krumholz et al.l[2005cl : 
iBonnell k Bate! 120061: lKrumholzll2006af) . which appears 
to be ruled out by observations of th e star formation 
rate and age spread iii young clusters (|Tan et al.l 120061 : 
iKrumholz fc Tanll2006D . The latter possibility requires 
stellar densities that would be very difficult to achieve 
without global collapse, and which are orders of magni- 
tude larger than the highest observed stellar densities in 
young clusters. 

Our results suggest that the solution to this dilemma 
is that massive stars do form directly from the collapse 
of massive cores (MT03), and that these cores do not 
fragment strongly because radiation feedback effectively 
shuts off fragmentation (|Krumholdi2006b( ) . In a massive, 
collapsing core, most of the mass goes into one primary 
object. Thus, massive cores - objects that we know exist 
from observations - are the direct progenitors of massive 
stars, with no need for intermediate steps of competitive 
accretion or collisions. 

In this work we have not addressed the question of 
whether, at higher masses, radiation pressure might halt 
accretion and thereby prevent massive cores from making 
massive stars. However, theoretical work to date suggests 
that radiation beaming by a combination of protostellar 
outflow cavit ies, the protostellar disk and the accret- 
ing en velope (jYorke fc Sonnhalterll2002l : IKrumholz et al.l 
I2005air0 ) provide a robust mechanism for allowing accre- 
tion to continue in the face of radiation pressure. We plan 
to report on three dimensional radiation-hydrodynamic 
simulations of this problem in future work. 

4.4. Implications for Future Simulations 

Our results show definitively that radiative transfer 
significantly modifies the manner in which accretion and 
fragmentation occur in the environments where massive 
stars form, at least on the size scales of individual cores. 
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The effective heating radius is > 1000 AU even before 
nuclear burning starts in any protostar, simply due to 
accretion luminosity. Once nuclear burning starts, the 
heating radius will rise rapidly, since the luminosity rises 
as roughly . It is not possible to capture this effect 
simply by using a modified equation of state, because the 
heating process depends on the radiative transfer of en- 
ergy from gas falling onto protostellar surfaces to gas in 
the surrounding envelope. We conclude that simulations 
of massive star formation and star cluster formation that 
do not include radiative transfer and accretion luminos- 
ity are not reliable on size scales below several thousand 
AU, and that such calculations almost certainly overes- 
timate the amount of fragmentation that occurs. This is 
true even before any of the stars formed begin nuclear 
burning. If one wishes to continue a simulation to the 
point where deuterium burning starts, one must include 
that effect as well. 

A critical question, which our work raises but does 
not address, is the extent to which the overfragmenta- 
tion problem affects simulations of low-mass star forma- 
tion. In such environments cores are usually separated 
by more than 1000 AU, accretion rates and the resulting 
accretion luminosities are lower, and lower column den- 
sities produce lower opacities to what radiation there is. 
Thus, we expect the effect on fragmentation to be less 
severe than we have found. Nonetheless, it seems likely 
that there will be some effect, particularly for models in 
which brown dw arfs or low-mass stars form by disk frag- 
mentation (e.g. iBate et al.l 12001 120031 : iGoodwin etall 
[2004) , and for competitive accretion models in which nu- 
merous brown dwarfs or low-mass stars form in clusters 
^ 1000 AU in size, and then evol ve in a manner dictated 
largely by N-body in teractions (iBate fc Bonnelll 120051: 
iBonnell fc Bate! |2005[ ). In particular. iMatzner fc LevinI 
(120051 ^ argue that radiative transfer effects are likely to 
prevent the formation of brown dwarfs by disk fragmen- 
tation, and our results support the idea that this effect 
might be important. It is therefore critical to repeat 
these calculations with radiative transfer to see if the 
fragments pe rsist once better phys i cs is i ncluded, a point 
also made by IWhitehouse fc Batd (|2006[ ). 

5. CONCLUSION 

We report the results of simulations of the collapse and 
fragmentation of massive protostellar cores using gravito- 
radiation hydrodynamics with adaptive mesh refinement. 
We find that including radiative transfer in our simula- 
tions produces dramatic effects on the evolution of these 
objects. When radiation is included, massive protostel- 
lar cores with the properties of observed cores do not 
fragment strongly. They collapse to a handful of ob- 
jects, with the majority of the mass accreting onto one 
primary object. The object gains mass by accretion of 
gas that is prevented from collapsing by radiative heat- 
ing. Some low-mass stars do form in addition to the 
primary massive star, through a combination of frag- 
mentation at sites sufficiently far from the protostar to 
be fairly cool, and fragmentation of unstable protostel- 



lar disks around the massive star. However, these do 
not contain significant mass. The disks that form are 
able to transport mass onto the central star very rapidly 
due to a large-scale gravitational instability, which ap- 
pears to form due to the SLING mechanism. Overall, 
our results are consistent with the turbulent core model 
of MT02 and MT03, with the analyti c treatment ofradia- 
tive suppression of fragmentation bv iKrumhol j (|2006bD . 
a nd with the analytic massive protostellar disk models 
of iKratter fc Matzneil ()2006[ ). 

Our results suggest that massive cores are the direct 
progenitors of massive stars, without an intermediate 
phase of competitive accretion or stellar collisions. Both 
of these mechanisms require that a large number of proto- 
stars form out of dense gas in clusters ^ 1000 AU in size 
or smaller. However, such strong fragmentation seems 
not to occur in the environments where massive stars 
form, because rapid accretion rapidly raises the gas tem- 
perature and prevents nucleation of small protostars. In- 
stead, most gas in a massive cores accretes onto the first 
object to form, and this object prevents other, compara- 
ble mass stars from forming after it. 

One additional conclusion from our work is that one 
must include radiative transfer in simulations in order to 
obtain the correct behavior on small scales. Prescriptions 
for the equation of state based on the barotropic approx- 
imation or on optically thin heating and cooling models 
severely underestimate gas temperatures because they do 
not, and cannot, include heating by accretion luminos- 
ity onto a central protostar. The high densities found in 
massive protostellar cores mean that this luminosity can 

reach thousands of Lq even for ^ 1 Mq protostars, an 
effect that cannot be neglected. 
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